# load in libraries
import equadratures as eq
import numpy as np
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import scipy.integrate as integrate

pressure assumption: $\frac{p}{p_0} = y = R^{-3}f_1$ ^[1]\ density assumption: $\frac{\rho}{\rho_0} = \psi$ ^[2]\ velocity assumption: $u = R^{-3/2}\phi_1$ ^[3]\ where $f_1$, $\psi$, and $\phi_1$ are functions of $\eta = \frac{r}{R}$ where $r$ is the radial position and $R$ is the radius of the shock wave

Using the Navier-Stokes equation: $\rho \frac{Du}{Dt} = -\nabla p + \rho g + \mu \nabla^2 u$ under the assumption of no gravitational or viscous effects we have $\rho \frac{Du}{Dt} = -\nabla p$, using the defintion of a materal derivative and the equations for denisty, pressure, and velocity from above we get ^[4]:

$$ \large \rho \left(\frac{\partial u}{\partial t} + (u \cdot \nabla)u\right) = -\frac{\partial p}{\partial r} \longrightarrow \rho \left(\frac{\partial u}{\partial t} + (u \cdot \nabla)u\right) = -\frac{p_0}{\rho}\frac{\partial y}{\partial r} $$$$ \large \left(\frac{\partial (R^{-3/2}\phi_1)}{\partial t} + R^{-3/2}\phi_1 \left(\frac{\partial (R^{-3/2} \phi_1)}{\partial r}\right)\right) = -\frac{p_0}{\rho_0\psi} \frac{\partial ( R^{-3} f_1)}{\partial r} $$$$ \large R^{-3/2}\frac{\partial \phi_1}{\partial t} - \frac{3}{2}R^{-5/2}\frac{dR}{dt}\phi_1 + R^{-3/2}\phi_1 \left(R^{-3/2}\frac{\partial \phi_1}{\partial r}\right) = -\frac{p_0}{\rho_0 \psi} R^{-3} \frac{\partial f_1}{\partial r} $$

Since $f_1$ and $\phi_1$ are functions of $\eta$ their partial derivatives can be expanded:

$$ \large R^{-3/2}\frac{\partial \phi_1}{\partial \eta}\frac{\partial \eta}{\partial t}- \frac{3}{2}R^{-5/2}\frac{dR}{dt}\phi_1 + R^{-3/2}\phi_1 \left(R^{-3/2}\frac{\partial \phi_1}{\partial \eta}\frac{\partial \eta}{\partial r}\right) = -\frac{p_0}{\rho_0 \psi} R^{-3} \frac{\partial f_1}{\partial \eta}\frac{\partial \eta}{\partial r} $$

Writing $\frac{\partial \phi_1}{\partial \eta}$ as $\phi_1'$ and $\frac{\partial f_1}{\partial \eta}$ as $f_1'$ and using $\frac{\partial \eta}{\partial t} = \frac{\partial r/R}{\partial t} = -\frac{r}{R^2} \frac{dR}{dt}$ and $\frac{\partial \eta}{\partial r}=\frac{\partial r/R}{\partial t}=\frac{1}{R}$ we get that:

$$ \large -R^{-3/2}\phi_1'\frac{r}{R^2}\frac{dR}{dt}- \frac{3}{2}R^{-5/2}\frac{dR}{dt}\phi_1 + R^{-3/2}\phi_1 \left(R^{-3/2}\phi_1'\frac{1}{R}\right) = -\frac{p_0}{\rho_0 \psi} R^{-3} f_1'\frac{1}{R} $$

Simplifying this we have ^[5]:

$$ \large -\left(\frac{3}{2}\phi_1 + \eta\phi_1'\right)R^{-5/2}\frac{dR}{dt} + R^{-4}\left(\phi_1\phi_1' + \frac{p_0}{\rho_0 \psi}f_1'\right) = 0 $$

Taylor claims in equation 6 that this can be satisfied by $\frac{dR}{dt}=AR^{-3/2}$ ^[6] where A is some constant and in the end we could simplify ^[5] to get ^[7]): However plugging that in would get:

$$ \large -\left(\frac{3}{2}\phi_1 + \eta\phi_1'\right)R^{-5/2}AR^{-3/2} + R^{-4}\left(\phi_1\phi_1' + \frac{p_0}{\rho_0 \psi}f_1'\right) = 0 $$$$ \large \longrightarrow -A\left(\frac{3}{2}\phi_1 + \eta\phi_1'\right)R^{-15/4} + R^{-4}\left(\phi_1\phi_1' + \frac{p_0}{\rho_0 \psi}f_1'\right) = 0 $$

Which would then not cancel out the factor of $R$

Continuing on and using the equation of continuity $\frac{\partial \rho}{\partial t} + u \frac{\partial \rho}{\partial r} + \rho\left(\frac{\partial u}{\partial r} + \frac{2u}{r}\right)=0$ ^[8] we can use the same equations for denisty and velocity:

$$ \large \rho_0 \frac{\partial \psi}{\partial t} + \rho_0R^{-3/2}\phi_1 \frac{\partial \psi}{\partial r} + \rho_0\psi \left(R^{-3/2}\frac{\partial \phi_1}{\partial r} + \frac{2R^{-3/2}\phi_1}{r}\right)=0 $$$$ \large \rho_0 \frac{\partial \psi}{\partial \eta}\frac{\partial \eta}{\partial t} + \rho_0R^{-3/2}\phi_1 \frac{\partial \psi}{\partial \eta}\frac{\partial \eta}{\partial r} + \rho_0\psi \left(R^{-3/2}\frac{\partial \phi_1}{\partial \eta}{\partial \eta}{\partial r} + \frac{2R^{-3/2}\phi_1}{r}\right)=0 $$$$ \large \rho_0 \psi'\left(-\frac{r}{R^2}\frac{dR}{dt}\right)+ \rho_0R^{-3/2}\phi_1 \psi'\frac{1}{R} + \rho_0\psi \left(R^{-3/2}\phi_1'\frac{1}{R} + \frac{2R^{-3/2}\phi_1}{r}\right)=0 $$$$ \large \rho_0 \psi'\left(-\frac{r}{R^2}AR^{-3/2}\right)+ \rho_0R^{-5/2}\phi_1 \psi' + \rho_0\psi \left(R^{-5/2}\phi_1' + \frac{2R^{-3/2}\phi_1}{r}\right)=0 $$$$ \large -A\eta \psi'R^{-5/2}+ R^{-5/2}\phi_1 \psi' + \psi \left(R^{-5/2}\phi_1' + \frac{2R^{-5/2}\phi_1}{\eta}\right)=0 $$

Simplifying this we have ^[9]:

$$ \large -A\eta \psi'+ \phi_1 \psi' + \psi \left(\phi_1' + \frac{2\phi_1}{\eta}\right)=0 \ (9) $$

For a perfect gas $p\rho^{-\gamma} = const$ therefore $\frac{D(p\rho^{-\gamma})}{Dt}=0$ and therefore $\frac{\partial (p\rho^{-\gamma})}{\partial t} + u \frac{\partial (p \rho^{-\gamma})}{\partial r} =0$ ^[10]

$$ \large p \frac{\partial \rho^{-\gamma}}{\partial t} + \rho^{-\gamma} \frac{\partial p}{\partial t} + u\left(p \frac{\partial \rho^{-\gamma}}{\partial r} + \rho^{-\gamma} \frac{\partial p}{\partial r}\right)=0 $$$$ \large -\gamma p \rho^{-\gamma-1}\frac{\partial \rho}{\partial t} + \rho^{-\gamma} \frac{\partial p}{\partial t} + u\left(-\gamma p \rho^{-\gamma-1}\frac{\partial \rho}{\partial r} + \rho^{-\gamma} \frac{\partial p}{\partial r}\right)=0 $$$$ \large -\gamma p \rho^{-1}\frac{\partial \rho}{\partial t} + \frac{\partial p}{\partial t} + u\left(-\gamma p \rho^{-1}\frac{\partial \rho}{\partial r} + \frac{\partial p}{\partial r}\right)=0 $$$$ \large -\gamma p_0R^{-3}f_1 \rho_0(\rho_0 \psi)^{-1}\frac{\partial \psi}{\partial t} -3 f_1 p_0 R^{-4}\frac{dR}{dt} + p_0R^{-3}\frac{\partial f_1}{\partial t} + u\left(-\gamma p_0R^{-3}f_1 \rho_0(\rho_0 \psi)^{-1}\frac{\partial \psi}{\partial r} + p_0R^{-3}\frac{\partial f_1}{\partial r}\right)=0 $$$$ \large -\gamma p_0R^{-3}f_1 \psi^{-1} \psi' \left(-\frac{r}{R^2}\frac{dR}{dt}\right) -3 f_1 p_0 R^{-4}\frac{dR}{dt} + R^{-3}p_0f_1' \left(-\frac{r}{R^2}\frac{dR}{dt}\right) + u\left(-\gamma p_0R^{-3}f_1 \psi^{-1}\psi'\frac{1}{R} + p_0R^{-3}f_1'\frac{1}{R}\right)=0 $$$$ \large A\gamma p_0R^{-11/2}f_1\eta \psi^{-1} \psi' -3 f_1 p_0 R^{-4}AR^{-3/2} + A\eta R^{-11/2}p_0f_1' + R^{-3/2}\phi_1 \left(-\gamma p_0R^{-4}f_1 \psi^{-1}\psi' + p_0R^{-4}f_1'\right)=0 $$$$ \large A\gamma f_1\eta \psi^{-1} \psi' -3f_1A + A\eta f_1' + \phi_1 \left(-\gamma f_1 \psi^{-1}\psi' + f_1'\right)=0 $$$$ \large -A(3f_1 + \eta f_1') +\gamma f_1\psi^{-1} \psi'( A\eta -\phi_1) + \phi_1f_1'=0 $$

Simplifying this becomes ^[11]: $$ \large A(3f_1 + \eta f_1') +\gamma f_1\frac{\psi'}{\psi}(\phi_1 - A\eta) - \phi_1f_1'=0 $$

We can reduce ^[7], ^[9], and ^[11] by forming nondimensional forms of $f_1$ and $\phi_1$ where $f=f_1\frac{a^2}{A^2}$ (12) and $\phi = \frac{\phi_1}{A}$ (13) where $a^2 = \gamma \frac{p_0}{\rho_0}$ and is the speed of sound. Doing this for ^[7] we get ^[7a]:

$$ \large -A\left(\frac{3}{2}A\phi + \eta A\phi'\right) + \left(A^2\phi\phi' + \frac{p_0}{\rho_0 \psi}\frac{A^2}{a^2}f'\right) = 0 $$$$ \large -\frac{3}{2}\phi + \frac{p_0}{\rho_0 \psi}\frac{\rho_0}{\gamma p_0}f' = \eta \phi' - \phi\phi' $$$$ \large \phi'(\eta - \phi)= \frac{1}{\gamma}\frac{f'}{\psi} -\frac{3}{2}\phi $$

Using the subsitutions in ^[9] gives ^[9a]:

$$ \large -A\eta \psi'+ A\phi \psi' + \psi \left(A\phi' + \frac{2A\phi}{\eta}\right)=0 $$$$ \large \psi \left(\phi' + \frac{2\phi}{\eta}\right)=\psi'(\eta -\phi) $$$$ \large \frac{\psi'}{\psi} = \frac{\left(\phi' + \frac{2\phi}{\eta}\right)}{\eta -\phi} $$

Using the substitution in ^[11] gives ^[11a]:

$$ \large A(3\frac{A^2}{a^2}f + \eta \frac{A^2}{a^2}f') +\gamma \frac{A^2}{a^2}f\frac{\psi'}{\psi}(A\phi - A\eta) - A\phi\frac{A^2}{a^2}f'=0 $$$$ \large 3f + \eta f' +\gamma f\frac{\psi'}{\psi}(\phi - \eta) - \phi f'=0 $$

Substituting ^[7a] and ^[9a] into ^[11a] gives ^[14]:

$$ \large 3f + \eta f' +\gamma f\frac{\left(\phi' + \frac{2\phi}{\eta}\right)}{\eta -\phi}(\phi - \eta) - \phi f'=0 $$$$ \large 3f + \eta f' -\gamma f\left(\phi' + \frac{2\phi}{\eta}\right) - \phi f'=0 $$$$ \large 3f + \eta f' -\gamma f\left(\frac{\frac{1}{\gamma}\frac{f'}{\psi} - \frac{3}{2}\phi}{\eta - \phi} + \frac{2\phi}{\eta}\right) - \phi f'=0 $$$$ \large f'\left((\eta - \phi)^2 - \frac{f}{\psi}\right) = f\left(-3\eta + \phi\left(3 + \frac{1}{2}\gamma\right) - \frac{2\gamma \phi^2}{\eta}\right) $$

Using the Rankine-Hugoniot relations and reducing we get ^[15], ^[16], and ^[17]:

$$ \large \frac{\rho_1}{\rho_0} = \frac{\gamma - 1 +(\gamma+1)y_1}{\gamma+1 + (\gamma-1)y_1} $$$$ \large \frac{U^2}{a^2} = \frac{1}{2\gamma}(\gamma -1 + (\gamma+1)y_1) $$$$ \large \frac{u_1}{U} = \frac{2(y_1-1)}{\gamma-1 + (\gamma+1)y_1} $$

Where $\rho_1$, $y_1$, and $u_1$ are the values of $\rho$, $\gamma$, and $y_1$ immediately behind/before the shock and $U=\frac{dR}{dt}$, we can simplify these for the case where $y_1 >>$ to get ^[15a], ^[16a], and ^[17a]:

$$ \large \frac{\rho_1}{\rho_0} = \frac{\gamma+1}{\gamma-1} $$$$ \large \frac{U^2}{a^2} = \frac{\gamma +1}{2\gamma}y_1 $$$$ \large \frac{u_1}{U} = \frac{2}{\gamma+1} $$

Since $\psi=\frac{\rho}{\rho_0}$ ^[15a] becomes ^[15b]:

$$ \large \psi = \frac{\gamma+1}{\gamma-1} $$

Next, since $f=\frac{a^2}{A^2}$ then $a^2=\frac{A^2f}{f_1}$ and because $U=\frac{dR}{dt}=AR^{-3/2}$ then $U^2=A^2R^{-3}$, additionally $y=R^{-3}f_1$ using these equations ^[16a] becomes ^[16b]:

$$ \large \frac{A^2R^{-3}f_1}{A^2f} = \frac{\gamma + 1}{2\gamma}R^{-3}f_1 $$$$ \large f = \frac{2\gamma}{\gamma + 1} $$

Lastly, for ^[17a] we can use that $u=R^{-3/2}A\phi$ and $U=\frac{dR}{dt}=AR^{-3/2}$ to get ^[17b]:

$$ \large \frac{AR^{-3/2}\phi}{AR^{-3/2}} = \frac{2}{\gamma+1} $$$$ \large \phi = \frac{2}{\gamma+1} $$

Note: In Taylor's paper ^[16a] is different: However, this does not derive from equation 16 and does not simplify to ^[16b] so I'm not sure what he did but the end result of ^[16b] is the same

By using (^[15b], ^[16b], ^[17b]) and using that $\gamma=1.4$ we can get values for $f$, $\phi$, and $\psi$ at $\eta=1$:

# define values at eta=1
f_0 = 7/6
phi_0 = 5/6
psi_0 = 6
g = 1.4

Using ^[14], ^[7a], and ^[9a] we can define functions to calculation $f'$, $\phi'$, and $\psi'$

# define f'
def dev_f(f, phi, psi, n):
    f_prime = (f*(-3*n+phi*(3+0.5*g)-(2*g*phi**2)/n))/((n-phi)**2 -f/psi)
    return f_prime 

# define phi'
def dev_phi(f, phi, psi, n):
    f_prime = dev_f(f, phi, psi, n)
    phi_prime = ((1/g)*(f_prime/psi)-(3/2)*(phi))/(n-phi)
    return phi_prime

# define psi'
def dev_psi(f, phi, psi, n):
    phi_prime = dev_phi(f, phi, psi, n)
    psi_prime = psi*(phi_prime+2*phi/n)/(n-phi)
    return psi_prime

To find an approximate function for $f$, $\phi$, and $\psi$ we can calcuate the change at very small increments

ns = np.linspace(1,0.5, 100)

# step by step calculation 
fs = []
phis = []
psis = []

f = f_0
phi = phi_0
psi = psi_0

f_prime_list = []

delta = np.mean(np.diff(ns))
for n in ns:
    fs.append(f)
    phis.append(phi)
    psis.append(psi)
    
    f_prime = dev_f(f, phi, psi, n)
    phi_prime = dev_phi(f, phi, psi, n)
    psi_prime = dev_psi(f, phi, psi, n) 
    f_prime_list.append(f_prime)
    
    f += delta*f_prime
    phi += delta*phi_prime
    psi += delta*psi_prime

The output for $\phi$ is very nearly linear and can be approximated using ^[14] and ^[9a]:

$$ \large 3f + \eta f' -\gamma f\left(\phi' + \frac{2\phi}{\eta}\right) - \phi f'=0 $$

Note: in Taylor ^[19] is written differently but both forms are equivalent \ This can be simplifed to get ^[20]

$$ \large \frac{f'}{f}(\eta - \phi) = \gamma \phi' + \frac{2\gamma \phi}{\eta} - 3 $$

If the left-hand side is ignored then the solution is $\phi = \frac{\eta}{\gamma}$ ^[21] this linear relationship follows what is calculated by the step process however it starts to break away at high values of $\eta$ meaning there is a correction factor of the form $\alpha\eta^n$ where $\alpha$ and n are constants ^[22]

$$ \large \phi = \frac{\eta}{\gamma} + \alpha\eta^n $$

These constants can be found using ^[16b] and ^[17b], firstly by using ^[17b] and plugging that value in for $\phi$ at $\eta=1$ we get ^[23]:

$$ \large \frac{2}{\gamma +1} = \frac{1}{\gamma} + \alpha (1)^n \longrightarrow \alpha = \frac{\gamma -1}{\gamma(\gamma +1)} $$

Next, we can use (^[15b], ^[16b], ^[17b]) to find the true value of $f'$ at $\eta=1$ getting that $f' = \frac{2\gamma^2 +7\gamma-3}{\gamma-1}$, next we can take the derivative of ^[22] to get $\phi' = \frac{1}{\gamma} +\alpha n \eta^{n-1}$ then plug this value into ^[20] and set them equal in order to calculate n ^[24]:

$$ \large n = \frac{7\gamma-1}{\gamma^2 -1} $$

Then we can use ^[22] to solve for $f$ and $\psi$, firstly plugging ^[22] into ^[20] we get ^[25]:

$$ \large \frac{f'}{f} = \frac{(n+2)\alpha\gamma^2\eta^{n-2}}{\gamma-1-\gamma\alpha\eta^{n-1}} $$

Taking the derivative and solving for the derivation constant using the value of $f$ at $\eta=1$ we get that ^[26]:

$$ \large ln(f) = ln\left(\frac{2\gamma}{\gamma+1}\right) - \frac{2\gamma^2+7\gamma-3}{7-\gamma}ln\left(\frac{\gamma+1-\eta^{n-1})}{\gamma}\right) $$

Next we can use ^[22] and its derivative to solve for $\psi$ by using ^[9a]:

$$ \large \frac{\psi'}{\psi} = \frac{\frac{1}{\gamma} + \alpha n \eta^{n-1} + \frac{2}{\eta}\left(\frac{\eta}{\gamma} + \alpha \eta^n\right)}{\eta - \left(\frac{\eta}{\gamma} + \alpha\eta^n\right)}=\frac{3+(n+2)\gamma\alpha\eta^{n-1}}{(\gamma-1)\eta-\alpha\gamma\eta^n} $$

Taking the derivative and solving for $\psi$ gets ^[28]:

$$ \large ln(\psi) = ln\left(\frac{\gamma+1}{\gamma-1}\right)+\frac{3}{\gamma-1}ln(\eta)-2\frac{\gamma+5}{7-\gamma}ln\left(\frac{\gamma+1-\eta^{n-1}}{\gamma}\right) $$

Below we can check these equations for $\phi$, $f$, and $\psi$

# guess constants and values
a = (g-1)/(g*(g+1))
N = (7*g-1)/(g**2-1)
phi_guess = ns/g
phi_guess_corrected = ns/g + a*ns**N
f_guess_corrected = (2*g/(g+1))*(((g+1)/g)-((ns**(N-1))/g))**(-(2*g**2+7*g-3)/(7-g))
psi_guess_corrected = ((g+1)/(g-1))*(ns**(3/(g-1)))*((g+1-ns**(N-1))/g)**(-2*(g+5)/(7-g))

# taylors data
taylor_ns = np.linspace(1,0.5,26)
taylor_fs = [f_0, 0.949, 0.808, 0.711, 0.643, 0.593, 0.556, 0.528, 0.507, 0.491, 0.478, 0.468, 0.461, 0.455, 0.45, 0.447, 
             0.444, 0.442, 0.44, 0.439, 0.438, 0.438, 0.437, 0.437, 0.437, 0.436]
taylor_phis = [phi_0, 0.798, 0.767, 0.737, 0.711, 0.687, 0.665, 0.644, 0.625, 0.607, 0.590, 0.573, 0.557, 0.542, 0.527, 
               0.513, 0.498, 0.484, 0.47, 0.456, 0.443, 0.428, 0.415, 0.402, 0.389, 0.375]
taylor_psis = [psi_0, 4, 2.808, 2.052, 1.534, 1.177, 0.919, 0.727, 0.578, 0.462, 0.37, 0.297, 0.239, 0.191, 0.152, 0.12, 
               0.095, 0.074, 0.058, 0.044, 0.034, 0.026, 0.019, 0.014, 0.01, 0.007]

fig = go.Figure()

fig.add_scatter(x=ns, y=fs, name='step', line_color = 'skyblue')
fig.add_scatter(x=taylor_ns, y=taylor_fs, name='Taylor', line_color='darkviolet', line_dash='dash')
fig.add_scatter(x=ns, y=f_guess_corrected, name='corrected guess', line_dash='dot', line_color='deeppink')

fig.add_scatter(x=ns, y=phis, visible=False, name='step', line_color='skyblue')
fig.add_scatter(x=taylor_ns, y=taylor_phis, visible=False, name='Taylor', line_dash='dash', line_color='darkviolet')
fig.add_scatter(x=ns, y=phi_guess, visible=False, name='guess', line_dash='dashdot', line_color='pink')
fig.add_scatter(x=ns, y=phi_guess_corrected, visible=False, name='corrected guess', line_dash='dot', line_color='deeppink')

fig.add_scatter(x=ns, y=psis, visible=False, name='step', line_color='skyblue')
fig.add_scatter(x=taylor_ns, y=taylor_psis, visible=False, name='Taylor', line_dash='dash', line_color='darkviolet')
fig.add_scatter(x=ns, y=psi_guess_corrected, visible=False, name='corrected guess', line_dash='dot', line_color='deeppink')

fig.update_layout(updatemenus=[dict(active=0, buttons=[dict(label='f', method='update', 
                  args=[{'visible':[True, True, True, False, False, False, False, False, False, False]}]), 
                  dict(label=r'$\phi$', method='update', args=[{'visible':[False, False, False, True, True, True, True, 
                  False, False, False]}]), dict(label=r'$\psi$', method='update', args=[{'visible':[False, False, False, 
                  False, False, False, False, True, True, True]}])])])

fig.update_xaxes(title=r'$\eta \left(\frac{r}{R}\right)$')

fig.show()

To calculate the energy of the explosion we need to look at two parts, kinetic and heat \ Kinetic energy in integral form using spherical coordinates is $KE=4\pi\int_0^Rfrac{1}{2}\rho u^2r^2dr$ \ In terms of the variables $\phi$, $\psi$, $f$, and $\eta$ we can write this as:

$$ \large 2\pi\int_0^11 \rho_0\psi A^2R^{-3}\phi^2 R^2\eta^2 Rd\eta = 2\pi A^2\rho_0\int_0^1 \psi \phi^2 \eta^2 d\eta $$

Heat energy in integral form using spherical coordinates is $HE=4\pi\int_0^R\frac{pr^2}{\gamma-1}dr$ \ In terms of the variables $\phi$, $\psi$, $f$, and $\eta$ we can write this as:

$$ \large \frac{4\pi}{\gamma-1} \int_0^1 p_0 \frac{A^2}{a^2}fR^{-3} R^2\eta^2 Rd\eta = \frac{4\pi A^2 \rho_0}{\gamma(\gamma-1)} \int_0^1 f \eta^2 d\eta $$

Therefore the total Energy of the system can be written as:

$$ \large E = 2\pi A^2 \rho_0 \left(\int_0^1 \psi \phi^2 \eta^2 d\eta + \frac{2}{\gamma(\gamma-1)}\int_0^1 f \eta^2 d\eta\right) $$

Using that $\frac{dR}{dt}=AR^{-3/2}$ we can integrate and solve for A to get that $A=\frac{2}{5}R^{5/2}t^{-1}$ therefore:

$$ \large E = 2\pi \frac{4}{25}R^5t^{-2} \rho_0 \left(\int_0^1 \psi \phi^2 \eta^2 d\eta + \frac{2}{\gamma(\gamma-1)}\int_0^1 f \eta^2 d\eta\right) $$

We can state that $K=\frac{4}{5}\left(2\pi\int_0^1 \psi \phi^2 \eta^2 d\eta + \frac{4\pi}{\gamma(\gamma-1)}\int_0^1 f \eta^2 d\eta\right)$ so that the energy is $E=K\rho_0R^5t^{-2}$

# define a function that calculates K based on gamma
def K(gamma):
    a = (gamma-1)/(gamma*(gamma+1))
    N = (7*gamma-1)/(gamma**2-1)
    D = ((gamma+1)/(gamma-1))*((gamma+1)/gamma)**(-2*(gamma+5)/(7-gamma))
    I_1 = (lambda x: (x**2)*((x/gamma + a*x**N)**2)*(D*x**(3/(gamma-1))))
    I_2 = (lambda x: (x**2)*(2*gamma/(gamma+1))*((gamma+1)/gamma - (x**(N-1))/gamma)**(-(2*gamma**2+7*gamma-3)/(7-gamma)))
    K = (4/25)*(2*np.pi*integrate.quad(I_1, 0, 1)[0] + 4*np.pi*(1/(gamma*(gamma-1)))*integrate.quad(I_2, 0, 1)[0])
    return K

gammas = np.linspace(1.2, 1.7, 100)

Ks = []
for gamma in gammas:
    Ks.append(K(gamma))

In both part I and II I did not get the same values for the first integral (kinetic energy) so there is some minor errors/differences in Taylor's paper at $\gamma=1.4$ he got this value: I am not getting that value instead it is around $0.087$ the same goes for the value of that integral at all other values of gamma, I'm using the approximate guess formulas that he created and evaluated it using an online calculator as well as scipy so I am unsure what I am doing wrong

rho_0 = 1.25 # kg/m^3
ergs = 1e7 # ergs/J
tnt = 4.25e16 # ergs/tons TNT

From that $\frac{dR}{dt}=AR^{-3/2}$ the relationship between R and t can be found by taking the integral: $t=\frac{2}{5}AR^{5/2}$ meaning that the time is proportional to some constant times $R^{5/2}$ with the values taken from the images of the explosion we can find this constant:

R = np.array([11.1, 19.9, 25.4, 28.8, 31.9, 34.2, 36.3, 38.9, 41, 42.8, 44.4, 46, 46.9, 48.7, 59, 61.1, 62.9, 64.3, 65.6,
             67.3, 106.5, 130, 145, 175, 185]) #m
t = np.array([0.1, 0.24, 0.38, 0.52, 0.66, 0.8, 0.94, 1.06, 1.22, 1.36, 1.5, 1.65, 1.79, 1.93, 3.26, 3.53, 3.8, 4.07, 4.34,
              4.61, 15, 25, 34, 53, 62])*(1e-3) #millisec -> sec

# find the approximate value of R^5 t^-2
R5t2 = np.mean((R**5)*(t**(-2)))

# find the total energy based on gamma and convert to tons TNT
E = np.array(Ks)*rho_0*R5t2
E_ergs = E*ergs
E_TNT = E_ergs / tnt

# plot the energy vs the values that Taylor got 
fig = go.Figure()

fig.add_scatter(x=gammas, y=E)
fig.add_scatter(x=gammas, y=E_TNT, yaxis='y2', line_color='black', name='Me')
fig.add_scatter(x=[1.2, 1.3, 1.4, 1.667], y=[34000, 22900, 16800, 9500], yaxis='y2', mode='markers', marker_color='blue',
                name='Taylor')

fig.update_layout(yaxis=dict(title='Energy (ergs)'), yaxis2=dict(side='right', title='Tons TNT'), showlegend=False,
                  xaxis=dict(range=[1.15,1.75]))

fig.show()

First test released an approximate amount of 20,000 tons of TNT (source: Wikipedia. ) )

Taylor states that the approximate temperature behind the shock wave correlates to a gamma value of 1.3, at 1.3 on the graph the value is ~20,000 tons TNT

The pressure ratio is $\frac{p}{p_0}=R^{-3}f_1=\frac{A^2}{a^2}R^{-3}f$ using that $a^2=\gamma\frac{p_0}{\rho_1}$ we can get that $p=p_0A^2\frac{\rho_0}{\gamma p_0}R^{-3}f$, next using that $A=\frac{2}{5}R^{5/2}t^{-1}$ the pressure is $p=\frac{\rho_0}{\gamma}\frac{4}{25}R^{5}t^{-2}R^{-3}f$ to find the max pressure at each value of gamma we need the highest value of f which happens at $\eta=1$ and is $f=\frac{2\gamma}{\gamma+1}$ so the pressure becomes $p_{max}=\frac{2\rho_0}{\gamma+1}\frac{4}{25}R^{5}t^{-2}R^{-3}$ \ The energy is equivalent to $E=K\rho_0 R^{5}t^{-2}$ so the pressure is $p_{max}=\frac{8}{25}\frac{E}{K(\gamma+1)}R^-3$ where $\frac{E}{K}$ is a constant value

p_z = []
for gamma in gammas:
    g_p = []
    for r in R:
        p_max = (rho_0*R5t2)*(8/25)*(1/(1+gamma))*r**(-3)
        g_p.append(p_max)
    p_z.append(g_p)

# plot the pressure vs. the radius of the shock and the gamma value
fig = go.Figure()

fig.add_contour(z=np.log10(p_z), x=R, y=gammas, contours_coloring='heatmap', line_width=0, colorscale='viridis',
                colorbar=dict(title='Max Pressure [log(P)]', titleside='right'))
fig.update_layout(xaxis=dict(title='R [m]'), yaxis=dict(title=r'$\gamma$'))

fig.show()

Energy and pressure evaluations are taken from Taylor Part I equations 34-38 and Taylor Part II equations 1-8 as well as tables 1-3